home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / FOUR1.DEM < prev    next >
Text File  |  1991-04-29  |  3KB  |  98 lines
  1. PROGRAM d12r1(input,output);
  2. (* driver for routine FOUR1 *)
  3. CONST
  4.    nn=32;
  5.    nn2=64; (* 2*nn *)
  6. TYPE
  7.    gldarray = ARRAY [1..nn2] OF real;
  8. VAR
  9.    ii,i,isign : integer;
  10.    data,dcmp : gldarray;
  11.  
  12. PROCEDURE prntft(data : gldarray; nn : integer);
  13. VAR
  14.    ii,mm,n : integer;
  15. BEGIN
  16.    writeln('n':4,'real(n)':13,'imag.(n)':13,'real(N-n)':12,'imag.(N-n)':13);
  17.    writeln (0:4,data[1]:14:6,data[2]:12:6,data[1]:12:6,data[2]:12:6);
  18.    mm := nn DIV 2;
  19.    FOR ii := 1 to mm DO BEGIN
  20.       n := 2*ii+1;
  21.       writeln (((n-1) DIV 2):4,data[n]:14:6,data[n+1]:12:6,
  22.          data[2*nn+2-n]:12:6,data[2*nn+3-n]:12:6)
  23.    END;
  24.    writeln (' press return to continue ...');
  25.    readln
  26. END;
  27.  
  28. (*$I MODFILE.PAS *)
  29. (*$I FOUR1.PAS *)
  30.  
  31. BEGIN
  32.    writeln('h(t) := real-valued even-function');
  33.    writeln('h(n) := h(N-n) and real?');
  34.    FOR ii := 1 to nn DO BEGIN
  35.       i := 2*ii-1;
  36.       data[i] := 1.0/(sqr((i-nn-1.0)/nn)+1.0);
  37.       data[i+1] := 0.0
  38.    END;
  39.    isign := 1;
  40.    four1(data,nn,isign);
  41.    prntft(data,nn);
  42.    writeln('h(t) := imaginary-valued even-function');
  43.    writeln('h(n) := h(N-n) and imaginary?');
  44.    FOR ii := 1 to nn DO BEGIN
  45.       i := 2*ii-1;
  46.       data[i+1] := 1.0/(sqr((i-nn-1.0)/nn)+1.0);
  47.       data[i] := 0.0
  48.    END;
  49.    isign := 1;
  50.    four1(data,nn,isign);
  51.    prntft(data,nn);
  52.    writeln('h(t) := real-valued odd-function');
  53.    writeln('h(n) := -h(N-n) and imaginary?');
  54.    FOR ii := 1 to nn DO BEGIN
  55.       i := 2*ii-1;
  56.       data[i] := ((i-nn-1.0)/nn)/(sqr((i-nn-1.0)/nn)+1.0);
  57.       data[i+1] := 0.0
  58.    END;
  59.    data[1] := 0.0;
  60.    isign := 1;
  61.    four1(data,nn,isign);
  62.    prntft(data,nn);
  63.    writeln('h(t) := imaginary-valued odd-function');
  64.    writeln('h(n) := -h(N-n) and real?');
  65.    FOR ii := 1 to nn DO BEGIN
  66.       i := 2*ii-1;
  67.       data[i+1] := ((i-nn-1.0)/nn)/(sqr((i-nn-1.0)/nn)+1.0);
  68.       data[i] := 0.0
  69.    END;
  70.    data[2] := 0.0;
  71.    isign := 1;
  72.    four1(data,nn,isign);
  73.    prntft(data,nn);
  74. (* transform, inverse-transform test *)
  75.    FOR ii := 1 to nn DO BEGIN
  76.       i := 2*ii-1;
  77.       data[i] := 1.0/(sqr(0.5*(i-nn-1.0)/nn)+1.0);
  78.       dcmp[i] := data[i];
  79.       data[i+1] := (0.25*(i-nn-1.0)/nn)
  80.             *exp(-sqr(0.5*(i-nn-1.0)/nn));
  81.       dcmp[i+1] := data[i+1]
  82.    END;
  83.    isign := 1;
  84.    four1(data,nn,isign);
  85.    isign := -1;
  86.    four1(data,nn,isign);
  87.    writeln;
  88.    writeln('double fourier transform:':33,'original data:':23);
  89.    writeln;
  90.    writeln('k':4,'real h(k)':14,'imag h(k)':13,
  91.       'real h(k)':17,'imag h(k)':12);
  92.    FOR ii := 1 to (nn DIV 2) DO BEGIN
  93.       i := 2*ii-1;
  94.       writeln(((i+1) DIV 2):4,dcmp[i]:14:6,dcmp[i+1]:12:6,
  95.          data[i]/nn:17:6,data[i+1]/nn:12:6)
  96.    END
  97. END.
  98.